Fracture Toughness of Metallic Glasses: Ductile-to-Brittle Transition? 
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Quantitative understanding of the fracture toughness of metallic glasses, including the associated 
ductile-to-brittle transitions, is not yet available. Here we use a simple model of plastic deformation 
in glasses, coupled to an advanced Eulerian level set formulation for solving complex free boundary 
problems, to calculate the fracture toughness of metallic glasses as a function of the degree of 
structural relaxation corresponding to different annealing times near the glass temperature. Our 
main result indicates the existence of an elasto-plastic crack tip instability for sufficiently relaxed 
glasses, resulting in a marked drop in the toughness, which we interpret as a ductile-to-brittle 
transition similar to experimental observations. 



The mechanical properties of glassy materials still pose 
challenges of great scientific and technological impor- 
tance. One such fundamental property is the fracture 
toughness - the ability of a material to resist failure in 
the presence of a crack [1 . Theoretically predicting the 
fracture toughness of materials, which is lacking in gen- 
eral, is a particularly pressing problem in the context 
of metallic glasses. Metallic glasses constitute a promis- 
ing new class of materials, possessing superior properties, 
whose usage in structural applications is severely limited 
by their relatively low fracture toughness [2H2]. 

Recent observations demonstrated a marked drop in 
the fracture toughness of metallic glasses as a function 
of composition and degree of structural relaxation (con- 
trolled through annealing near the glass temperature T g ) 
[6l I9HTT] . The drop in the toughness, which is commonly 
correlated with Poisson's ratio [9UTT]. is interpreted as a 
kind of ductile-to-brittle transition [9j [11] . To the best 
of our knowledge, no basic theoretical understanding of 
this important observation is currently available. 

In this Letter we calculate the fracture toughness of 
metallic glasses based on the low-temperature Shear- 
Transformation-Zone (STZ) model, using an advanced 
Eulerian level set formulation for solving complex free 
boundary problems. We demonstrate the existence of 
an elasto-plastic crack tip instability as a function of in- 
creasing degree of structural relaxation, which results in 
a drop in the fracture toughness. We propose to interpret 
this instability as a ductile-to-brittle transition similar to 
the one observed experimentally. 

The Shear- Transformation- Zone (STZ) model of amor- 
phous plasticity [T2UT5] has recently been shown 
to emerge within a systematic formulation of non- 
equilibrium thermodynamics [I6j [17] and to capture a 
wide range of glassy deformation phenomena [15, 18H22]. 
Its main advantage in the present context is that it of- 
fers a way to quantify the degree of structural relaxation 
and the deformation-driven evolution of structural dis- 
order. Our goal here is to use the STZ model in a way 
that goes beyond previous analyses; rather than fixing 
the model parameters to quantitatively describe a given 



phenomenon, we treat it as a predictive model where its 
parameters are estimated from independent sources and 
another phenomenon - crack initiation - is studied. 

We focus here on a simple version of the STZ model, 
retaining only salient physical ingredients. As we are 
interested in the fracture toughness of metallic glasses 
at temperatures well below the glass temperature T g , we 
neglect all spontaneous, non-driven, relaxation processes 
and set the plastic rate of deformation D pl to zero for 
stresses below the shear yield stress s y . For s > s y we 
have 

DP l (s,T, X ) = To 1 A(x)C(s,T) [1 - s y /s] s/s , (1) 

where s = a— |trcr 1 is the deviatoric stress tensor (cr is 
the Cauchy stress) and \/2s = y/sjjsij [14] [15]. D pl is 
expressed as a product of physically meaningful terms. 
Tq -1 is a molecular vibration rate. A(x) is the probability 
to find a structural fluctuation that is particularly sus- 
ceptible to shear-driven rearrangements - an STZ. It is 
a function of an effective disorder temperature x> to be 
discussed below. C(s, T) quantifies the (dimensionless) 
rate in which STZs actually undergo shear transforma- 
tions as a function of stress and temperature. The last 
terms represent deformation-induced anisotropy ("back 
stress") and also make the whole expression tensorially 
consistent. 

The effective temperature x characterizes the out-of- 
equilibrium structural degrees of freedom of a glass [16] . 
It satisfies an effective heat equation of the form p~3j [16] 

roX = r(s,x) (Xoo -x) , (2) 

where again spontaneous thermally activated relaxation 
is excluded. Xoo is the steady state value of x and T(s, x) 
is a dimensionless strength of mechanically-generated 
noise that tends to rejuvenate the glass (when X<Xoo)- 
The very same theoretical framework predicts that A(x) 
m Eq. Q is given by a generalized Boltzmann factor, 
A(x) = ex P ( — e z/^BX)i where e z is a typical STZ for- 
mation energy. It is this physical description of struc- 
tural disorder that makes the STZ model most suitable 



2 



for studying the fracture toughness as a function of the 
degree of structural relaxation. 

To complete the presentation of the model in Eqs. 
Q-Q, we need to specify explicit forms for C(s, T) 
and r(s, x). The latter has been proposed to be pro- 
portional to the rate of plastic work DfjjSij [15 , i.e. 

T(s,x)=r Df j s ij /s y . C(s,T) = ±[K(s,T)+K(-s,T)] is 
the average of forward and backward STZ transition rates 
ft(±s,T)=exp (- A ^ T e °^ ), which we assume to follow 
a linearly stress-biased thermal activation process. Here 
A is the typical energy activation barrier, Q is the typi- 
cal activation volume and eo is the typical local strain at 
the transition [6j [23] . In the presence of the high stresses 
near a tip of a crack, Qe^s may become larger than A, in 
which case we assume the exponential thermal activation 
form crosses over to a much weaker dependence associ- 
ated with a linear, non-activated, dissipative mechanism 
[201 . Hence, 



(3) 



. f e~ A / kBT cosh [ft e s/k B T] forfte 5<A 
J [Qe s/2A for Qe s > A) 

As A> ksT, the two expressions connect continuously 
(but not differentiably ) . The slope of the linear relation 
was chosen so as not to introduce additional parameters. 
These details do not affect the qualitative nature of the 
results to follow. 

To proceed, we adopt an Eulerian formulation and 
write the total rate of deformation tensor as a sum of 
elastic and plastic contributions, D tot = D elj cD p \ where 
D tot = \[Vv+(Vv) T ), D el =d t e J r v • Ve+e • u—u • e and 
u) = ^[Vv — (Vt?) T ]. The strain tensor e is related to a 
through Hooke's law cr = K trel + 2p (e — |trel), where 
K and p are the bulk and shear moduli, respectively. 
The velocity field v(r,t), where r is the spatial coordi- 
nate, evolves through the momentum balance equation 
po (ptv+v- Vv) — V-cr, where po is the mass density (as- 
sumed constant hereafter). 

Consider a straight notch (crack) with root radius p 
(see Fig. [I]) under plane-strain conditions. A polar co- 
ordinate system (r, 6) is set a distance p/b behind the 
notch root and = is the symmetry axis. We adopt a 
boundary layer formulation in which the following univer- 
sal mode I (tensile) crack tip velocity fields are imposed 
on a scale much larger than p [24H29] 



v x {r,6,t) 
v y (r,0,t) 
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where Ki(t) is the mode I stress intensity factor and v is 
Poisson's ratio pQ. The main advantage of this approach 
is that the stress intensity factor uniquely couples the 
inner scales near the tip to the outer scales and hence can 



be controlled independently without solving the global 
crack problem pQ. 

The linear elastic fracture toughness is the critical 
value of the stress intensity factor, Kj ci at which the crack 
initiates and global failure occurs. There is ample exper- 
imental and numerical evidence that metallic glasses un- 
der tension fail locally near crack tips by the nucleation 
of voids [30 - 34 . We interpret this at the continuum level 
(atomistic aspects might be also relevant [35]) as a local 
cavitation instability initiating at a structural fluctuation 
when the hydrostatic tension |trcr exceeds a threshold, 
which for non-hardening materials is estimated as [36] 



2s y (l + log [2£/(3V3s y )]) /V3 , 



(5) 



where E is Young's modulus and s y /E<^l. 

The model parameters for Vitreloy 1, a widely studied 
metallic glass for which the annealing time dependence 
of the fracture toughness was measured [lOj [11] , are es- 
timated from independent sources, we set p = 37GPa, 
v = 0.35, po ^ 6g/cm 3 and s y ~ 0.85GPa [1 [6]. The 
basic vibrational timescale is To — 10 _13 s. The activa- 
tion volume of an STZ was estimated to be Q ~ 1000A 
[37] and typically e ~ 0.1 [6 , hence ^e ~ 100 A 3 . The 
typical activation barrier is of the order of leV; we set 
A = 0.7eV [37] . The STZ formation energy should be 
somewhat larger than A and we choose e z — 1.8eV. Fi- 
nally, the steady state value of the effective temperature 
is expected to be between T g = 623K and the melting tem- 
perature T m ~ 1000K. Previous works suggest Xoo — 900K 
[20]. We set T = 400K, well below T g . 

We set K / = 10MPav / ms- 1 and p = 65pm [10] E]. A 
key parameter is the initial value of the effective tem- 
perature, x( r ,t = Q)=Xo- In [10] E]: Vitreloy 1 was an- 
nealed for different times at T g and the fracture toughness 
dropped by an order of magnitude, from Ki c ~ 85MP^L^/m 
for the as-cast samples to Kj c ~ S.SMPa^/m for the 12 
hours annealed samples. Within the model, we represent 
the effect of increasing annealing times by decreasing val- 
ues of the initial effective temperature, and focus on the 
range xo = 600 — 660K. All other parameters remained 
fixed. 

We numerically solved the equations for cr, v and \ us ~ 
ing the recently proposed Eulerian finite-difference sim- 
ulation framework, where free boundaries are implicitly 
tracked by the level set method [38]. A key advantage 
of this method is its ability to naturally handle topolog- 
ical changes, such as those involved in material failure. 
The combination of finite-difference and level set meth- 
ods provides a flexible platform to study complex physi- 
cal phenomena such as crack initiation and propagation. 

The widely separated timescales of elastic and plas- 
tic deformations make our equations stiff. In [22] . an 
explicit update procedure, in which the timestep was 
chosen to be small enough to resolve elastic waves, was 
employed. It would be prohibitively computationally ex- 



3 




FIG. 1. (Color online) The normalized hydrostatic pres- 
sure field p/sy near a notch at various loading levels Ki — 
20, 40, 6OMPa0n (from top to bottom) for xo = 600K (panels 
a-c) and xo = 660K (panels d-f). A movie is available at [43] . 



pensive to access physically relevant timescales using this 
procedure. We therefore constructed a new numerical 
scheme in which a and x are explicitly updated, but v 
is solved for implicitly using quasi-static force balance 
V-cr = 0. Details of the quasi-static scheme and its veri- 
fication will be given elsewhere. Here we just stress that 
this scheme allows us to use physically realistic loading 
rates and to dynamically switch to the explicit scheme 
when rapid failure initiates. The calculations presented 
here employed a —20 < xj p, yj p < 20 domain, using a 
1025x 1025 grid. Increasing grid resolution and/or do- 
main size did not significantly affect the results. 

In Fig. [I] we plot a sequence of three snapshots of the 
hydrostatic pressure field p(r, t) = — §tnr for xo = 600K 
(more relaxed) and xo = 660K (less relaxed), taken at 
the same value of Ki. The two sequences seem to ex- 
hibit a similar qualitative behavior in which p attains 
a minimum ahead of the notch root at a distance that 
increases with Kj [26 29 . There are, however, marked 



quantitative differences; the lower xo exhibits a signifi- 
cantly smaller minimum (accompanied by a sharp spatial 
variation) and the local notch root radius of curvature 
decreases, suggesting the onset of a localization process. 

To further explore the crack tip dynamics, we plot in 
Fig. [5] two snapshots of the effective temperature x( r ^) 
for each xo- Recall that x( r ^) quantifies structural dis- 
order - the higher x> the higher the disorder and the 
easier it is to flow. Both the spatial distribution of x 
and the notch geometry are markedly different in the two 
cases. In the higher xo case, x( r ^) is rather smoothly 
distributed in the near tip region and the notch under- 
goes continuous blunting - its radius of curvature grows 
continuously and uniformly with Kj. 

The lower xo case is qualitatively different. Initially, 
at small loads, there is little plastic deformation and x 
remains nearly constant at its initial value xo- As Kj 
increases, plastic deformation localizes in the notch root 
vicinity, resulting in a more sharply and inhomogeneously 
distributed x? featuring small scale filamentary struc- 
tures. These dynamics are strongly coupled to the notch 
geometry; the radius of curvature of the notch varies spa- 
tially, with a pronounced reduction near the root. It is 
this localization process - an elasto-plastic crack tip in- 
stability - that is responsible for the marked differences 
in the minima of p in Fig. [I] 
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FIG. 2. (Color online) The effective temperature field x a t 
Kt 40MPaVm (top) and Kt SOMPa^m (bottom) for 
Xo=600K (panels a-b) and xo = 660K (panels c-d). A movie 
is available at [43] . 

What are the implications of this instability for the 
fracture toughness? As discussed above, large \p\ may 
induce the nucleation of voids, which might lead to catas- 
trophic failure. Therefore, we focus our attention on the 
minimum of the pressure Pmi n (t) =Min{p(r, t)}. In Fig. 
[3k we plot Pmin vs. Kj for the two Xo' s - At small Kj 
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both samples respond linear elastically (and hence identi- 
cally). As Kj increases, local near tip yielding occurs and 
the curves progressively and significantly deviate from 
the elastic line. Already here we observe some quantita- 
tive differences: the lower xo sample exhibits less plastic 
deformation and consequently less stress relaxation and 
tip blunting, resulting in more negative Pmin- As Kj 
further increases, a clear signature of the tip instability 
discussed above is observed, where p m i n drops abruptly 
for the lower xo? while the curve for the higher one ex- 
hibits smooth and moderate variation with Kj. 
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FIG. 3. (Color online) (a) Pmin /s y vs. Ki for xo = 600K (solid 
blue line) and xo = 660K (dashed red line). The horizontal 
line at p m in/s y = —4.5 is the threshold for void nucleation. 
(inset) Xmax vs. Ki (b) A snapshot of the system (p/s y is 
plotted) when a void (small white circle) nucleates, (c) The 
subsequent catastrophic failure. A movie is available at [43] . 

A complementary view on the elasto-plastic nature 
of the instability is obtained by plotting Xmax(f) = 
Max{x(f% £)}, which quantifies the magnitude of plas- 
tic deformation, vs. Kj in the inset of Fig. [3^i. For 
Xo = 660K, a linear elastic regime (Xmax — Xo) is fol- 
lowed by a smooth and moderate increase of Xmax to- 
ward Xoo- For xo = 600K the linear elastic regime is 
followed by an accelerated and very sharp increase of 
Xmax, which is mirrored in the drop of p m in in the main 
panel. To make things quantitative, we use E / s y 85 



for Vitreloy 1 in Eq. (5j to get a c ~ 5^ [39 . In Fig. 
[3^i, we chose 4.5s y as the threshold (horizontal line) for 
void nucleation, which suggests a large difference in the 
fracture toughness, Kj c ~ SOMPa^/m for xo = 600K and 
Kj c ~ 80MPay / m for xo = 660K. Varying <j c will not 
change the qualitative nature of this main result, though 
the flatness of the xo = 660K curve suggests quantitative 
implications. 

Does the nucleation of a void lead to catastrophic fail- 
ure? i.e. can we interpret Kj at which p m i n meets the 
threshold as the fracture toughness Ki c ? To address this 
issue we take advantage of the model's dynamical nature 
and the numerical method's flexibility to study the post 
void nucleation dynamics. A void nucleation is shown in 
Fig. [3J3- The subsequent dynamics, a snapshot of which 
is shown in panel (c), proceed through a rapid succession 
of void nucleations, leading to the coalescence of the ini- 
tial void with the notch root and to rapid crack propaga- 
tion which results in catastrophic failure. The emerging 
crack pattern is reminiscent of some experimental obser- 
vations (random fluctuations in the void nucleation loca- 
tions were introduced to avoid artificial grid effects) [40] . 
In light of this catastrophic failure, we interpret the large 
variation in K\ at which the threshold is met in Fig. [3^i 
for the two different xo's as a kind of ductile-to-brittle 
transition similar to the experimental observations. 

The crack tip instability, which leads to the marked 
drop in the fracture toughness discussed above, has both 
constitutive and geometric origins. The central physi- 
cal question here is how efficiently a material can tame 
the linear elastic stress singularity, associated with the 
universal crack tip fields of Eqs. Q, by stress relaxation 
processes [41 . Stress relaxation is mediated both by bulk 
plastic deformation and by the accompanying geometri- 
cal changes in the shape of the notch - the higher the ra- 
dius of curvature, the lower the stress concentration. As 
a glass becomes progressively more structurally relaxed 
(less disordered), these stress relaxation processes be- 
come progressively more limited and below some thresh- 
old a tip instability sets in. 

As mentioned above, the ductile-to-brittle transition 
is commonly correlated with Poisson's ratio v [6, 9] [IT]. 
We suspect that this correlation might not be deep, but 
rather represents the fact that both the elastic and plas- 
tic responses of a glass depend on its state of disorder, 
quantified here by x- Hence, while there should exist a 
configurational equation of state of the form v(x) [HI HI] 
(which was neglected in our calculations), its effect on the 
fracture toughness is expected to be secondary compared 
to the strong exponential dependence of D pl on e z /kBXo 
through A(x). Indeed, for our parameters A drops by 
more than an order of magnitude when xo decreases from 
660K to 600K. 

The typical fracture toughness values that emerge from 
our calculations seem to be in the right ballpark, with- 
out fine-tuning the model's parameters. We would not, 
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however, take this to imply that the present model has 
quantitative predictive powers as there are still uncertain- 
ties about the details of the model (e.g. the form of the 
rate factor) and the values of the parameters. On the 
other hand, we do advocate the view that the model can 
be used to qualitatively predict new phenomena, such as 
the crack tip instability discussed above. 
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